## Code to create Figure 3a and Figure 3b
## Stability of Immigration Preferences
## Alexander Kustov, Dillon Laaker, and Cassidy Reller
remove(list = ls())

library("survey")
library("foreign")

location <- "/Users/dillonlaaker/Box/Stability/"
setwd(paste(location, "data/datasets", sep = ""))
ines <- read.dta("ines_data.dta")
bes <- read.dta("bes_data.dta")
taps <- read.dta("taps_data.dta")
ncp <- read.dta("ncp_data.dta")
shp <- read.dta("shp_data.dta")
liss <- read.dta("liss_data.dta")
gles <- read.dta("gles_data.dta")
cces <- read.dta("cces_data.dta")
vsg <- read.dta("vsg_data.dta")

##INES

means1 <- apply(ines, 2, mean, na.rm=TRUE)[6:9]
means2 <- apply(ines, 2, mean, na.rm=TRUE)[10:13]
survey <- replicate(4, "INES")
monthssince <- c(12, 24, 48, 60)
inesmeans<- data.frame(survey, monthssince, means1, means2)

means1 <- apply(ines, 2, function(x) mean(x[ines$all==1], na.rm=TRUE))[6:9]
means2 <- apply(ines, 2, function(x) mean(x[ines$all==1], na.rm=TRUE))[10:13]
means1 <- c(i.same23,i.same24,i.same26,i.same27)
means2 <- c(i2.same23,i2.same24,i2.same26,i2.same27)
survey <- replicate(4, "INES - All Waves")
monthssince <- c(12, 24, 48, 60)
inesmeansa<- data.frame(survey, monthssince, means1, means2)

##BES

means11 <- apply(bes, 2, mean, na.rm=TRUE)[80:86]
means21 <- apply(bes, 2, mean, na.rm=TRUE)[87:93]
means12 <- apply(bes, 2, mean, na.rm=TRUE)[94:100]
means22 <- apply(bes, 2, mean, na.rm=TRUE)[101:107]
means13 <- apply(bes, 2, mean, na.rm=TRUE)[108:114]
means23 <- apply(bes, 2, mean, na.rm=TRUE)[115:121]
survey <- replicate(7, "BES")
monthssince <- c(3, 7, 13, 26, 27, 33, 38)
besmeans <- data.frame(survey, monthssince, means11, means21, means12, means22, means13, means23)
besmeans$mean1 <- apply(besmeans[,c(3,5,7)], 1, function(x) mean(x, na.rm=TRUE))
besmeans$mean2 <- apply(besmeans[,c(4,6,8)], 1, function(x) mean(x, na.rm=TRUE))

means11 <- apply(bes, 2, function(x) mean(x[bes$all==1], na.rm=TRUE))[80:86]
means21 <- apply(bes, 2, function(x) mean(x[bes$all==1], na.rm=TRUE))[87:93]
means12 <- apply(bes, 2,function(x) mean(x[bes$all==1], na.rm=TRUE))[94:100]
means22 <- apply(bes, 2, function(x) mean(x[bes$all==1], na.rm=TRUE))[101:107]
means13 <- apply(bes, 2, function(x) mean(x[bes$all==1], na.rm=TRUE))[108:114]
means23 <- apply(bes, 2, function(x) mean(x[bes$all==1], na.rm=TRUE))[115:121]
survey <- replicate(7, "BES - All Waves")
besmeansa <- data.frame(survey, monthssince, means11, means21, means12, means22, means13, means23)
besmeansa$mean1 <- apply(besmeansa[,c(3,5,7)], 1, function(x) mean(x, na.rm=TRUE))
besmeansa$mean2 <- apply(besmeansa[,c(4,6,8)], 1, function(x) mean(x, na.rm=TRUE))

##TAPS

means1 <- apply(taps, 2, mean, na.rm=TRUE)[12:21]
means2 <- replicate(10, NA)
survey <- replicate(10, "TAPS")
monthssince <- c(7,14,19,26,29,38,41,43,46,48)
tapsmeans<- data.frame(survey, monthssince, means1, means2)

means1 <- apply(taps, 2, function(x) mean(x[taps$all==1], na.rm=TRUE))[12:21]
survey <- replicate(10, "TAPS - All Waves")
monthssince <- c(7,14,19,26,29,38,41,43,46,48)
tapsmeansa<- data.frame(survey, monthssince, means1, means2)

##LISS

means11 <- apply(liss, 2, mean, na.rm=TRUE)[180:187]
means21 <- apply(liss, 2, mean, na.rm=TRUE)[188:195]
means12 <- apply(liss, 2, mean, na.rm=TRUE)[196:203]
means22 <- apply(liss, 2, mean, na.rm=TRUE)[204:211]
means13 <- apply(liss, 2, mean, na.rm=TRUE)[212:219]
means23 <- apply(liss, 2, mean, na.rm=TRUE)[220:227]
means14 <- apply(liss, 2, mean, na.rm=TRUE)[228:235]
means24 <- apply(liss, 2, mean, na.rm=TRUE)[236:243]
means15 <- apply(liss, 2, mean, na.rm=TRUE)[244:251]
means25 <- apply(liss, 2, mean, na.rm=TRUE)[252:259]
means16 <- apply(liss, 2, mean, na.rm=TRUE)[260:267]
means26 <- apply(liss, 2, mean, na.rm=TRUE)[268:275]

survey <- replicate(8, "LISS")
monthssince <- c(12, 24, 36, 48, 60, 72, 96, 108)
lissmeans <- data.frame(survey, monthssince, means11, means21, means12, means22, means13, means23, means14, means24, means15, means25, means16, means26)
lissmeans$mean1 <- apply(lissmeans[,c(3,5,7,9,11,13)], 1, function(x) mean(x, na.rm=TRUE))
lissmeans$mean2 <- apply(lissmeans[,c(4,6,8,10,12,14)], 1, function(x) mean(x, na.rm=TRUE))

means11 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[180:187]
means21 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[188:195]
means12 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[196:203]
means22 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[204:211]
means13 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[212:219]
means23 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[220:227]
means14 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[228:235]
means24 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[236:243]
means15 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[244:251]
means25 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[252:259]
means16 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[260:267]
means26 <- apply(liss, 2, function(x) mean(x[liss$all==1], na.rm=TRUE))[268:275]

survey <- replicate(8, "LISS")
lissmeansa <- data.frame(survey, monthssince, means11, means21, means12, means22, means13, means23, means14, means24, means15, means25, means16, means26)
lissmeansa$mean1 <- apply(lissmeansa[,c(3,5,7,9,11,13)], 1, function(x) mean(x, na.rm=TRUE))
lissmeansa$mean2 <- apply(lissmeansa[,c(4,6,8,10,12,14)], 1, function(x) mean(x, na.rm=TRUE))


##ncp
n.sameimm1_34<-mean(ncp$sameimm1_34, na.rm=TRUE)
n.sameimm1_35<-mean(ncp$sameimm1_35, na.rm=TRUE)
n.sameimm1_36<-mean(ncp$sameimm1_36, na.rm=TRUE)
n.sameimm1_37<-mean(ncp$sameimm1_37, na.rm=TRUE)
n.sameimm1_38<-mean(ncp$sameimm1_38, na.rm=TRUE)
n2.sameimm1_34<-mean(ncp$same2imm1_34, na.rm=TRUE)
n2.sameimm1_35<-mean(ncp$same2imm1_35, na.rm=TRUE)
n2.sameimm1_36<-mean(ncp$same2imm1_36, na.rm=TRUE)
n2.sameimm1_37<-mean(ncp$same2imm1_37, na.rm=TRUE)
n2.sameimm1_38<-mean(ncp$same2imm1_38, na.rm=TRUE)
means11<-c(n.sameimm1_34, n.sameimm1_35, n.sameimm1_36, n.sameimm1_37, n.sameimm1_38)
means21<-c(n2.sameimm1_34, n2.sameimm1_35, n2.sameimm1_36, n2.sameimm1_37, n2.sameimm1_38)
survey <- replicate(5, "NCP - Question 1")
monthssince <- c(5, 13, 17, 25, 29)
ncp1means <- data.frame(survey, monthssince, means11, means21)

n.sameimm1_34<-mean(ncp$sameimm1_34[ncp$all==1], na.rm=TRUE)
n.sameimm1_35<-mean(ncp$sameimm1_35[ncp$all==1], na.rm=TRUE)
n.sameimm1_36<-mean(ncp$sameimm1_36[ncp$all==1], na.rm=TRUE)
n.sameimm1_37<-mean(ncp$sameimm1_37[ncp$all==1], na.rm=TRUE)
n.sameimm1_38<-mean(ncp$sameimm1_38[ncp$all==1], na.rm=TRUE)
n2.sameimm1_34<-mean(ncp$same2imm1_34[ncp$all==1], na.rm=TRUE)
n2.sameimm1_35<-mean(ncp$same2imm1_35[ncp$all==1], na.rm=TRUE)
n2.sameimm1_36<-mean(ncp$same2imm1_36[ncp$all==1], na.rm=TRUE)
n2.sameimm1_37<-mean(ncp$same2imm1_37[ncp$all==1], na.rm=TRUE)
n2.sameimm1_38<-mean(ncp$same2imm1_38[ncp$all==1], na.rm=TRUE)
means11<-c(n.sameimm1_34, n.sameimm1_35, n.sameimm1_36, n.sameimm1_37, n.sameimm1_38)
means21<-c(n2.sameimm1_34, n2.sameimm1_35, n2.sameimm1_36, n2.sameimm1_37, n2.sameimm1_38)
survey <- replicate(5, "NCP - Question 1 - All Waves")
monthssince <- c(5, 13, 17, 25, 29)
ncp1meansall <- data.frame(survey, monthssince, means11, means21)

n.sameimm2_34<-mean(ncp$sameimm2_34, na.rm=TRUE)
n.sameimm2_35<-mean(ncp$sameimm2_35, na.rm=TRUE)
n.sameimm2_36<-mean(ncp$sameimm2_36, na.rm=TRUE)
n.sameimm2_37<-mean(ncp$sameimm2_37, na.rm=TRUE)
n.sameimm2_38<-mean(ncp$sameimm2_38, na.rm=TRUE)
n2.sameimm2_34<-mean(ncp$same2imm2_34, na.rm=TRUE)
n2.sameimm2_35<-mean(ncp$same2imm2_35, na.rm=TRUE)
n2.sameimm2_36<-mean(ncp$same2imm2_36, na.rm=TRUE)
n2.sameimm2_37<-mean(ncp$same2imm2_37, na.rm=TRUE)
n2.sameimm2_38<-mean(ncp$same2imm2_38, na.rm=TRUE)
means12<-c(n.sameimm2_34, n.sameimm2_35, n.sameimm2_36, n.sameimm2_37, n.sameimm2_38)
means22<-c(n2.sameimm2_34, n2.sameimm2_35, n2.sameimm2_36, n2.sameimm2_37, n2.sameimm2_38)
survey <- replicate(5, "NCP - Question 2")
monthssince <- c(5, 13, 17, 25, 29)
ncp2means <- data.frame(survey, monthssince, means12, means22)

n.sameimm2_34<-mean(ncp$sameimm2_34[ncp$all==1], na.rm=TRUE)
n.sameimm2_35<-mean(ncp$sameimm2_35[ncp$all==1], na.rm=TRUE)
n.sameimm2_36<-mean(ncp$sameimm2_36[ncp$all==1], na.rm=TRUE)
n.sameimm2_37<-mean(ncp$sameimm2_37[ncp$all==1], na.rm=TRUE)
n.sameimm2_38<-mean(ncp$sameimm2_38[ncp$all==1], na.rm=TRUE)
n2.sameimm2_34<-mean(ncp$same2imm2_34[ncp$all==1], na.rm=TRUE)
n2.sameimm2_35<-mean(ncp$same2imm2_35[ncp$all==1], na.rm=TRUE)
n2.sameimm2_36<-mean(ncp$same2imm2_36[ncp$all==1], na.rm=TRUE)
n2.sameimm2_37<-mean(ncp$same2imm2_37[ncp$all==1], na.rm=TRUE)
n2.sameimm2_38<-mean(ncp$same2imm2_38[ncp$all==1], na.rm=TRUE)
means12<-c(n.sameimm2_34, n.sameimm2_35, n.sameimm2_36, n.sameimm2_37, n.sameimm2_38)
means22<-c(n2.sameimm2_34, n2.sameimm2_35, n2.sameimm2_36, n2.sameimm2_37, n2.sameimm2_38)
survey <- replicate(5, "NCP - Question 2 - All Waves")
monthssince <- c(7, 13, 17, 25, 29)
ncp2meansall <- data.frame(survey, monthssince, means12, means22)

ncpmeans1 <- data.frame(ncp1means, ncp2means)
ncpmeans1$means1 <- (ncpmeans1$means11 + ncpmeans1$means12)/2
ncpmeans1$means2 <- (ncpmeans1$means21 + ncpmeans1$means22)/2
survey <- replicate(5, "NCP")
monthssince <- c(7, 13, 17, 25, 29)
ncpmeans <- data.frame(survey, monthssince)
ncpmeans$means1 <- ncpmeans1$means1
ncpmeans$means2 <- ncpmeans1$means2

ncpmeans1 <- data.frame(ncp1meansall, ncp2meansall)
ncpmeans1$means1 <- (ncpmeans1$means11 + ncpmeans1$means12)/2
ncpmeans1$means2 <- (ncpmeans1$means21 + ncpmeans1$means22)/2
survey <- replicate(5, "NCP - All Waves")
monthssince <- c(7, 13, 17, 25, 29)
ncpmeansall <- data.frame(survey, monthssince)
ncpmeansall$means1 <- ncpmeans1$means1
ncpmeansall$means2 <- ncpmeans1$means2

##SHP
s.sameimm990<-mean(shp$sameimm990, na.rm=TRUE)
s.sameimm991<-mean(shp$sameimm991, na.rm=TRUE)
s.sameimm992<-mean(shp$sameimm992, na.rm=TRUE)
s.sameimm993<-mean(shp$sameimm993, na.rm=TRUE)
s.sameimm994<-mean(shp$sameimm994, na.rm=TRUE)
s.sameimm995<-mean(shp$sameimm995, na.rm=TRUE)
s.sameimm996<-mean(shp$sameimm996, na.rm=TRUE)
s.sameimm997<-mean(shp$sameimm997, na.rm=TRUE)
s.sameimm998<-mean(shp$sameimm998, na.rm=TRUE)
s.sameimm999<-mean(shp$sameimm999, na.rm=TRUE)
s.sameimm9911<-mean(shp$sameimm9911, na.rm=TRUE)
means1<-c(s.sameimm990, s.sameimm991, s.sameimm992, s.sameimm993, s.sameimm994, s.sameimm995, s.sameimm996, s.sameimm997, s.sameimm998, s.sameimm999, s.sameimm9911)
means2 <- replicate(11, NA)
survey <- replicate(11, "SHP")
monthssince <- c(12,24,36,48,60,72,84,96,108,120,144)
shpmeans <- data.frame(survey, monthssince, means1, means2)

s.sameimm990<-mean(shp$sameimm990[shp$all==1], na.rm=TRUE)
s.sameimm991<-mean(shp$sameimm991[shp$all==1], na.rm=TRUE)
s.sameimm992<-mean(shp$sameimm992[shp$all==1], na.rm=TRUE)
s.sameimm993<-mean(shp$sameimm993[shp$all==1], na.rm=TRUE)
s.sameimm994<-mean(shp$sameimm994[shp$all==1], na.rm=TRUE)
s.sameimm995<-mean(shp$sameimm995[shp$all==1], na.rm=TRUE)
s.sameimm996<-mean(shp$sameimm996[shp$all==1], na.rm=TRUE)
s.sameimm997<-mean(shp$sameimm997[shp$all==1], na.rm=TRUE)
s.sameimm998<-mean(shp$sameimm998[shp$all==1], na.rm=TRUE)
s.sameimm999<-mean(shp$sameimm999[shp$all==1], na.rm=TRUE)
s.sameimm9911<-mean(shp$sameimm9911[shp$all==1], na.rm=TRUE)
means1 <- c(s.sameimm990, s.sameimm991, s.sameimm992, s.sameimm993, s.sameimm994, s.sameimm995, s.sameimm996, s.sameimm997, s.sameimm998, s.sameimm999, s.sameimm9911)
means2 <- replicate(11, NA)
survey <- replicate(11, "SHP - All Waves")
monthssince <- c(12,24,36,48,60,72,84,96,108,120,144)
shpmeansall <- data.frame(survey, monthssince, means1, means2)

##GLES

g.sameimm12 <- mean(gles$sameimm12, na.rm=TRUE)
g.sameimm14 <- mean(gles$sameimm14, na.rm=TRUE)
g.sameimm16 <- mean(gles$sameimm16, na.rm=TRUE)
g.sameimm110 <- mean(gles$sameimm110, na.rm=TRUE)
g.sameimm111 <- mean(gles$sameimm111, na.rm=TRUE)
g.sameimm112 <- mean(gles$sameimm112, na.rm=TRUE)
g.sameimm113 <- mean(gles$sameimm113, na.rm=TRUE)
g.sameimm115 <- mean(gles$sameimm115, na.rm=TRUE)
g.sameimm117 <- mean(gles$sameimm117, na.rm=TRUE)
g.same2imm12 <- mean(gles$same2imm12, na.rm=TRUE)
g.same2imm14 <- mean(gles$same2imm14, na.rm=TRUE)
g.same2imm16 <- mean(gles$same2imm16, na.rm=TRUE)
g.same2imm110 <- mean(gles$same2imm110, na.rm=TRUE)
g.same2imm111 <- mean(gles$same2imm111, na.rm=TRUE)
g.same2imm112 <- mean(gles$same2imm112, na.rm=TRUE)
g.same2imm113 <- mean(gles$same2imm113, na.rm=TRUE)
g.same2imm115 <- mean(gles$same2imm115, na.rm=TRUE)
g.same2imm117 <- mean(gles$same2imm117, na.rm=TRUE)

means1 <- c(g.sameimm12, g.sameimm14, g.sameimm16, g.sameimm110, g.sameimm111, g.sameimm112, g.sameimm113, g.sameimm115, g.sameimm117)
means2 <- c(g.same2imm12, g.same2imm14, g.same2imm16, g.same2imm110, g.same2imm111, g.same2imm112, g.same2imm113, g.same2imm115, g.same2imm117)
survey <- replicate(9, "GLES")
monthssince <- c(1,2,3,40,45,48,50,52,53)
glesmeans <- data.frame(survey, monthssince, means1, means2)

g.sameimm12 <- mean(gles$sameimm12[gles$all==1], na.rm=TRUE)
g.sameimm14 <- mean(gles$sameimm14[gles$all==1], na.rm=TRUE)
g.sameimm16 <- mean(gles$sameimm16[gles$all==1], na.rm=TRUE)
g.sameimm110 <- mean(gles$sameimm110[gles$all==1], na.rm=TRUE)
g.sameimm111 <- mean(gles$sameimm111[gles$all==1], na.rm=TRUE)
g.sameimm112 <- mean(gles$sameimm112[gles$all==1], na.rm=TRUE)
g.sameimm113 <- mean(gles$sameimm113[gles$all==1], na.rm=TRUE)
g.sameimm115 <- mean(gles$sameimm115[gles$all==1], na.rm=TRUE)
g.sameimm117 <- mean(gles$sameimm117[gles$all==1], na.rm=TRUE)
g.same2imm12 <- mean(gles$same2imm12[gles$all==1], na.rm=TRUE)
g.same2imm14 <- mean(gles$same2imm14[gles$all==1], na.rm=TRUE)
g.same2imm16 <- mean(gles$same2imm16[gles$all==1], na.rm=TRUE)
g.same2imm110 <- mean(gles$same2imm110[gles$all==1], na.rm=TRUE)
g.same2imm111 <- mean(gles$same2imm111[gles$all==1], na.rm=TRUE)
g.same2imm112 <- mean(gles$same2imm112[gles$all==1], na.rm=TRUE)
g.same2imm113 <- mean(gles$same2imm113[gles$all==1], na.rm=TRUE)
g.same2imm115 <- mean(gles$same2imm115[gles$all==1], na.rm=TRUE)
g.same2imm117 <- mean(gles$same2imm117[gles$all==1], na.rm=TRUE)

means1 <- c(g.sameimm12, g.sameimm14, g.sameimm16, g.sameimm110, g.sameimm111, g.sameimm112, g.sameimm113, g.sameimm115, g.sameimm117)
means2 <- c(g.same2imm12, g.same2imm14, g.same2imm16, g.same2imm110, g.same2imm111, g.same2imm112, g.same2imm113, g.same2imm115, g.same2imm117)
survey <- replicate(9, "GLES - All Waves")
monthssince <- c(1,2,3,40,45,48,50,52,53)
glesmeansall <- data.frame(survey, monthssince, means1, means2)

##CCES

c.sameimm1_1012 <- mean(cces$sameimm1_1012, na.rm=TRUE)
c.sameimm1_1014 <- mean(cces$sameimm1_1014, na.rm=TRUE)
c.sameimm2_1012 <- mean(cces$sameimm2_1012, na.rm=TRUE)
c.sameimm2_1014 <- mean(cces$sameimm2_1014, na.rm=TRUE)
c.sameimm3_1012 <- mean(cces$sameimm3_1012, na.rm=TRUE)
c.sameimm3_1014 <- mean(cces$sameimm3_1014, na.rm=TRUE)
means11 <- c(c.sameimm1_1012, c.sameimm1_1014)
means21 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 1")
monthssince <- c(24, 48)
cces1means <- data.frame(survey, monthssince, means11, means21)
means12 <- c(c.sameimm2_1012, c.sameimm2_1014)
means22 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 2")
monthssince <- c(24, 48)
cces2means <- data.frame(survey, monthssince, means12, means22)
means13 <- c(c.sameimm3_1012, c.sameimm3_1014)
means23 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 3")
monthssince <- c(24, 48)
cces3means <- data.frame(survey, monthssince, means13, means23)

ccesmeans1 <- data.frame(cces1means, cces2means, cces3means)
ccesmeans1$means1 <- (ccesmeans1$means11 + ccesmeans1$means12 + ccesmeans1$means13)/3
ccesmeans1$means2 <- NA
survey <- replicate(2, "CCES")
monthssince <- c(24, 48)
ccesmeans <- data.frame(survey, monthssince)
ccesmeans$means1 <- ccesmeans1$means1
ccesmeans$means2 <- ccesmeans1$means2

c.sameimm1_1012 <- mean(cces$sameimm1_1012[cces$all==1], na.rm=TRUE)
c.sameimm1_1014 <- mean(cces$sameimm1_1014[cces$all==1], na.rm=TRUE)
c.sameimm2_1012 <- mean(cces$sameimm2_1012[cces$all==1], na.rm=TRUE)
c.sameimm2_1014 <- mean(cces$sameimm2_1014[cces$all==1], na.rm=TRUE)
c.sameimm3_1012 <- mean(cces$sameimm3_1012[cces$all==1], na.rm=TRUE)
c.sameimm3_1014 <- mean(cces$sameimm3_1014[cces$all==1], na.rm=TRUE)
means11 <- c(c.sameimm1_1012, c.sameimm1_1014)
means21 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 1 - All Waves")
monthssince <- c(24, 48)
cces1means <- data.frame(survey, monthssince, means11, means21)
means12 <- c(c.sameimm2_1012, c.sameimm2_1014)
means22 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 2 - All Waves")
monthssince <- c(24, 48)
cces2means <- data.frame(survey, monthssince, means12, means22)
means13 <- c(c.sameimm3_1012, c.sameimm3_1014)
means23 <- replicate(2, NA)
survey <- replicate(2, "CCES - Question 3 - All Waves")
monthssince <- c(24, 48)
cces3means <- data.frame(survey, monthssince, means13, means23)

ccesmeans1 <- data.frame(cces1means, cces2means, cces3means)
ccesmeans1$means1 <- (ccesmeans1$means11 + ccesmeans1$means12 + ccesmeans1$means13)/3
ccesmeans1$means2 <- NA
survey <- replicate(2, "CCES - All Waves")
monthssince <- c(24, 48)
ccesmeansall <- data.frame(survey, monthssince)
ccesmeansall$means1 <- ccesmeans1$means1
ccesmeansall$means2 <- ccesmeans1$means2

##vsg

c.sameimm1_1116 <- mean(vsg$sameimm1_1116, na.rm=TRUE)
c.sameimm1_1118 <- mean(vsg$sameimm1_1118, na.rm=TRUE)
c.sameimm2_1116 <- mean(vsg$sameimm2_1116, na.rm=TRUE)
c.sameimm2_1118 <- mean(vsg$sameimm2_1118, na.rm=TRUE)
c.sameimm3_1116 <- mean(vsg$sameimm3_1116, na.rm=TRUE)
c.sameimm3_1118 <- mean(vsg$sameimm3_1118, na.rm=TRUE)
means11 <- c(c.sameimm1_1116, c.sameimm1_1118)
means21 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 1")
monthssince <- c(60, 77)
vsg1means <- data.frame(survey, monthssince, means11, means21)
means12 <- c(c.sameimm2_1116, c.sameimm2_1118)
means22 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 2")
monthssince <- c(60, 77)
vsg2means <- data.frame(survey, monthssince, means12, means22)
means13 <- c(c.sameimm3_1116, c.sameimm3_1118)
means23 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 3")
monthssince <- c(60, 77)
vsg3means <- data.frame(survey, monthssince, means13, means23)

vsgmeans1 <- data.frame(vsg1means, vsg2means, vsg3means)
vsgmeans1$means1 <- (vsgmeans1$means11 + vsgmeans1$means12 + vsgmeans1$means13)/3
vsgmeans1$means2 <- NA
survey <- replicate(2, "VSG")
monthssince <- c(60, 77)
vsgmeans <- data.frame(survey, monthssince)
vsgmeans$means1 <- vsgmeans1$means1
vsgmeans$means2 <- vsgmeans1$means2

c.sameimm1_1116 <- mean(vsg$sameimm1_1116[vsg$all==1], na.rm=TRUE)
c.sameimm1_1118 <- mean(vsg$sameimm1_1118[vsg$all==1], na.rm=TRUE)
c.sameimm2_1116 <- mean(vsg$sameimm2_1116[vsg$all==1], na.rm=TRUE)
c.sameimm2_1118 <- mean(vsg$sameimm2_1118[vsg$all==1], na.rm=TRUE)
c.sameimm3_1116 <- mean(vsg$sameimm3_1116[vsg$all==1], na.rm=TRUE)
c.sameimm3_1118 <- mean(vsg$sameimm3_1118[vsg$all==1], na.rm=TRUE)
means11 <- c(c.sameimm1_1116, c.sameimm1_1118)
means21 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 1 - All Waves")
monthssince <- c(60, 77)
vsg1means <- data.frame(survey, monthssince, means11, means21)
means12 <- c(c.sameimm2_1116, c.sameimm2_1118)
means22 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 2 - All Waves")
monthssince <- c(60, 77)
vsg2means <- data.frame(survey, monthssince, means12, means22)
means13 <- c(c.sameimm3_1116, c.sameimm3_1118)
means23 <- replicate(2, NA)
survey <- replicate(2, "VSG - Question 3 - All Waves")
monthssince <- c(60, 77)
vsg3means <- data.frame(survey, monthssince, means13, means23)

vsgmeans1 <- data.frame(vsg1means, vsg2means, vsg3means)
vsgmeans1$means1 <- (vsgmeans1$means11 + vsgmeans1$means12 + vsgmeans1$means13)/3
vsgmeans1$means2 <- NA
survey <- replicate(2, "VSG - All Waves")
monthssince <- c(60, 77)
vsgmeansall <- data.frame(survey, monthssince)
vsgmeansall$means1 <- vsgmeans1$means1
vsgmeansall$means2 <- vsgmeans1$means2

figure3 <- rbind(inesmeans, besmeans, lissmeans, tapsmeans, ncpmeans, shpmeans, glesmeans, ccesmeans, vsgmeans, inesmeansall, besmeansall, lissmeansall, tapsmeansall, ncpmeansall, shpmeansall, glesmeansall, ccesmeansall, vsgmeansall)

setwd("/Users/dillonlaaker/Box Sync/Stability/Draft")
pdf("figure3a.pdf", height=6, width=12)
par(mgp=c(3,1,0))
par(mar=c(5,5,.75,1.5)) 
plot(figure3$monthssince[figure3$survey=="TAPS"], figure3$means1[figure3$survey=="TAPS"], type="o", lty=1, pch=4, cex=1.25, bg="white", ylim=c(0,1), yaxt="n", xaxt="n", bty="l", xaxs="i", yaxs="i", xlim=c(0,150), ylab="", xlab="")
mtext(side=2, text="Percentage giving same response", line=3, cex=1.25)
axis(side=2, at=c(0,.2,.4,.6,.8,1), lty=1, las=2)
axis(side=1, at=c(0, 25, 50, 75, 100, 125, 150), lty=1, las=0)
lines(figure3$monthssince[figure3$survey=="LISS"], figure3$means1[figure3$survey=="LISS"], type="o", lty=1, pch=22, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="SHP"], figure3$means1[figure3$survey=="SHP"], type="o", lty=1, pch=9, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="INES"], figure3$means1[figure3$survey=="INES"], type="o", lty=1, pch=24, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="BES"], figure3$means1[figure3$survey=="BES"], type="o", lty=1, pch=23, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="NCP"], figure3$means1[figure3$survey=="NCP"], type="o", lty=1, pch=21, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="GLES"], figure3$means1[figure3$survey=="GLES"], type="o", lty=1, pch=25, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="CCES"], figure3$means1[figure3$survey=="CCES"], type="o", lty=1, pch=11, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="VSG"], figure3$means1[figure3$survey=="VSG"], type="o", lty=1, pch=8, bg="white", cex=1)
legend("bottomright", legend=c("BES (7 cats.)", "CCES (2 cats.)", "GLES (7 cats.)", "INES (7 cats.)", "LISS (5 cats.)", "NCP (7 cats.)", "SHP (3 cats.)", "TAPS (2 cats.)", "VSG (3 cats.)"), pch=c(23, 11, 25, 24, 22, 21, 9, 4, 8), bty="n", pt.cex=1.25, cex=1.25)
dev.off()

pdf("figure3b.pdf", height=6, width=12) 
par(mgp=c(3,1,0))
par(mar=c(5,5,1,1.5))
plot(figure3$monthssince[figure3$survey=="GLES"], figure3$means2[figure3$survey=="GLES"], type="o", lty=1, pch=21, cex=1, bg="white", ylim=c(0,1), yaxt="n", xaxt="n", bty="l", xaxs="i", yaxs="i", xlim=c(0,150), ylab="", xlab="")
axis(side=2, at=c(0,.2,.4,.6,.8,1), lty=1, las=2)
axis(side=1, at=c(0, 25, 50, 75, 100, 125, 150), lty=1, las=0)
mtext(side=2, text="Percent changing by less than 2 categories", line=3, cex=1.25)
mtext(side=1, text="Months since first wave", line=3, cex=1.25)
lines(figure3$monthssince[figure3$survey=="LISS"], figure3$means2[figure3$survey=="LISS"], type="o", lty=1, pch=22, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="INES"], figure3$means2[figure3$survey=="INES"], type="o", lty=1, pch=24, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="BES"], figure3$means2[figure3$survey=="BES"], type="o", lty=1, pch=23, bg="white", cex=1)
lines(figure3$monthssince[figure3$survey=="NCP"], figure3$means2[figure3$survey=="NCP"], type="o", lty=1, pch=21, bg="white", cex=1)
legend("bottomright", legend=c("BES (7 cats.)", "GLES (7 cats.)", "INES (7 cats.)", "LISS (5 cats.)", "NCP (7 cats.)"), pch=c(23, 25, 24, 22, 21), bty="n", pt.cex=1.25, cex=1.25)
dev.off()